Load Data

dataset <- read.delim("raw_data/FigureS5N.txt", stringsAsFactors = FALSE)

dataset$genotype <- gsub(" ", "", dataset$genotype )
dataset$genotype <- factor(dataset$genotype)
dataset$Experiment <- factor(rep(paste0("exp", 1:(length(dataset$genotype)/length(levels(dataset$genotype)))), each=length(unique(dataset$genotype))))

dataset$UID <- factor(paste(dataset$Experiment, dataset$genotype))

# wide format
kable(dataset, row.names = F)
genotype NT H2O2_5uM H2O2_30uM H2O2_50uM Experiment UID
WT 1097 575 397 106 exp1 exp1 WT
YFP-ALC1 1425 1018 740 100 exp1 exp1 YFP-ALC1
WT 1768 942 462 95 exp2 exp2 WT
YFP-ALC1 1271 931 528 150 exp2 exp2 YFP-ALC1
WT 1252 741 456 95 exp3 exp3 WT
YFP-ALC1 1354 935 624 48 exp3 exp3 YFP-ALC1
library(reshape2)
# reshape to long format
dataset <- melt(dataset, variable.name = "Treatment", value.name = "Counts")

dataset$genotype <- relevel(dataset$genotype, ref = "WT")
dataset$UID <- relevel(dataset$UID, ref = "exp1 WT")

dataset$H2O2 <- gsub("NT","1",dataset$Treatment)
dataset$H2O2 <- gsub("H2O2_|uM","",dataset$H2O2)
dataset$H2O2 <- log10(as.integer(dataset$H2O2))




dataset$Offset <- NA
for(uid in levels(dataset$UID)){
        dataset$Offset[dataset$UID == uid] <- mean(dataset$Counts[dataset$UID == uid])
}

dataset$NormCounts <- dataset$Counts / dataset$Offset



dataset$Offset2 <- NA
for(gid in levels(dataset$genotype)){
        dataset$Offset2[dataset$genotype == gid] <- mean(dataset$NormCounts[dataset$genotype == gid & dataset$H2O2 == 0])
}

dataset$NormCounts2 <- dataset$NormCounts / dataset$Offset2



# long format
kable(dataset, row.names = F)
genotype Experiment UID Treatment Counts H2O2 Offset NormCounts Offset2 NormCounts2
WT exp1 exp1 WT NT 1097 0.000000 543.75 2.0174713 2.050234 0.9840200
YFP-ALC1 exp1 exp1 YFP-ALC1 NT 1425 0.000000 820.75 1.7362169 1.776869 0.9771216
WT exp2 exp2 WT NT 1768 0.000000 816.75 2.1646771 2.050234 1.0558196
YFP-ALC1 exp2 exp2 YFP-ALC1 NT 1271 0.000000 720.00 1.7652778 1.776869 0.9934767
WT exp3 exp3 WT NT 1252 0.000000 636.00 1.9685535 2.050234 0.9601604
YFP-ALC1 exp3 exp3 YFP-ALC1 NT 1354 0.000000 740.25 1.8291118 1.776869 1.0294017
WT exp1 exp1 WT H2O2_5uM 575 0.698970 543.75 1.0574713 2.050234 0.5157808
YFP-ALC1 exp1 exp1 YFP-ALC1 H2O2_5uM 1018 0.698970 820.75 1.2403290 1.776869 0.6980419
WT exp2 exp2 WT H2O2_5uM 942 0.698970 816.75 1.1533517 2.050234 0.5625464
YFP-ALC1 exp2 exp2 YFP-ALC1 H2O2_5uM 931 0.698970 720.00 1.2930556 1.776869 0.7277158
WT exp3 exp3 WT H2O2_5uM 741 0.698970 636.00 1.1650943 2.050234 0.5682739
YFP-ALC1 exp3 exp3 YFP-ALC1 H2O2_5uM 935 0.698970 740.25 1.2630868 1.776869 0.7108498
WT exp1 exp1 WT H2O2_30uM 397 1.477121 543.75 0.7301149 2.050234 0.3561130
YFP-ALC1 exp1 exp1 YFP-ALC1 H2O2_30uM 740 1.477121 820.75 0.9016144 1.776869 0.5074175
WT exp2 exp2 WT H2O2_30uM 462 1.477121 816.75 0.5656566 2.050234 0.2758985
YFP-ALC1 exp2 exp2 YFP-ALC1 H2O2_30uM 528 1.477121 720.00 0.7333333 1.776869 0.4127110
WT exp3 exp3 WT H2O2_30uM 456 1.477121 636.00 0.7169811 2.050234 0.3497070
YFP-ALC1 exp3 exp3 YFP-ALC1 H2O2_30uM 624 1.477121 740.25 0.8429585 1.776869 0.4744067
WT exp1 exp1 WT H2O2_50uM 106 1.698970 543.75 0.1949425 2.050234 0.0950831
YFP-ALC1 exp1 exp1 YFP-ALC1 H2O2_50uM 100 1.698970 820.75 0.1218398 1.776869 0.0685699
WT exp2 exp2 WT H2O2_50uM 95 1.698970 816.75 0.1163147 2.050234 0.0567324
YFP-ALC1 exp2 exp2 YFP-ALC1 H2O2_50uM 150 1.698970 720.00 0.2083333 1.776869 0.1172474
WT exp3 exp3 WT H2O2_50uM 95 1.698970 636.00 0.1493711 2.050234 0.0728556
YFP-ALC1 exp3 exp3 YFP-ALC1 H2O2_50uM 48 1.698970 740.25 0.0648430 1.776869 0.0364928

Plot Data

library(ggplot2)

# raw data
ggplot(dataset, aes(x=H2O2, y=Counts)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE, aes(colour=genotype)) +
        geom_point(aes(colour=genotype, shape=Experiment), size=2) +        
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_shape_manual(values=15:20) +
        scale_color_manual(values=c("#000000","#808000"))

# NormCounts Linear
ggplot(dataset, aes(x=H2O2, y=NormCounts, color=genotype)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=genotype), size=2) +        
        geom_smooth(method=lm, formula = y ~ x, se=FALSE) +
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_color_manual(values=c("#000000","#808000"))

# NormCounts2 Linear
ggplot(dataset, aes(x=H2O2, y=NormCounts2, color=genotype)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=genotype), size=2) +        
        geom_smooth(method=lm, formula = y ~ x, se=FALSE) +
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_color_manual(values=c("#000000","#808000"))

# NormCounts Quadratic
ggplot(dataset, aes(x=H2O2, y=NormCounts, color=genotype)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=genotype), size=2) +        
        geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE) +
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)")+
        scale_color_manual(values=c("#000000","#808000"))

# NormCounts2 Quadratic
ggplot(dataset, aes(x=H2O2, y=NormCounts2, color=genotype)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=genotype), size=2) +        
        geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE) +
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_color_manual(values=c("#000000","#808000"))

# NormCounts Cubic
ggplot(dataset, aes(x=H2O2, y=NormCounts, color=genotype)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=genotype), size=2) +        
        geom_smooth(method=lm, formula = y ~ poly(x,3), se=FALSE) +
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)")+
        scale_color_manual(values=c("#000000","#808000"))

# NormCounts2 Cubic
ggplot(dataset, aes(x=H2O2, y=NormCounts2, color=genotype)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=genotype), size=2) +        
        geom_smooth(method=lm, formula = y ~ poly(x,3), se=FALSE) +
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_color_manual(values=c("#000000","#808000"))

library(Cairo)

cairo_pdf("FigureS5N.pdf", width = 5, height = 4, family = "Arial")

ggplot(dataset, aes(x=H2O2, y=NormCounts2)) + 
        theme_bw() +
        theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(), 
              axis.line = element_line(colour = "black"), text = element_text(size=14),
              panel.border = element_blank(), panel.background = element_blank()) +
        geom_point(aes(colour = genotype)) +
        geom_smooth(method=lm, formula = y ~ poly(x,3), se=TRUE, aes(colour = genotype), fill='#CCCCCC') +
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        ylab(label = "Normalized Counts") +
        scale_color_manual(values=c("#000000","#808000"))

dev.off()
## quartz_off_screen 
##                 2

Models

library(MASS)
library(DHARMa)
library(lme4)
library(lmerTest)
library(bbmle)

Linear formula

fit1 <- lm(Counts ~ Experiment + H2O2*genotype, data = dataset)
print(summary(fit1))
## 
## Call:
## lm(formula = Counts ~ Experiment + H2O2 * genotype, data = dataset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -245.04 -132.06  -19.98   95.30  382.01 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           1299.866    110.649  11.748 7.10e-10 ***
## Experimentexp2          86.125     94.856   0.908    0.376    
## Experimentexp3           5.875     94.856   0.062    0.951    
## H2O2                  -686.474     81.573  -8.416 1.18e-07 ***
## genotypeYFP-ALC1        60.660    135.971   0.446    0.661    
## H2O2:genotypeYFP-ALC1   35.276    115.361   0.306    0.763    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 189.7 on 18 degrees of freedom
## Multiple R-squared:  0.8839, Adjusted R-squared:  0.8517 
## F-statistic: 27.42 on 5 and 18 DF,  p-value: 7.994e-08
cat("AIC: ", AIC(fit1))
## AIC:  326.9888
simres <- simulateResiduals(fittedModel = fit1)
plot(simres)

fit2 <- lm(NormCounts ~ H2O2*genotype, data = dataset)
print(summary(fit2))
## 
## Call:
## lm(formula = NormCounts ~ H2O2 * genotype, data = dataset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30870 -0.11268 -0.01016  0.10500  0.33774 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            1.98687    0.09414  21.105 3.86e-15 ***
## H2O2                  -1.01869    0.07987 -12.754 4.60e-11 ***
## genotypeYFP-ALC1      -0.15575    0.13314  -1.170    0.256    
## H2O2:genotypeYFP-ALC1  0.16077    0.11295   1.423    0.170    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1858 on 20 degrees of freedom
## Multiple R-squared:  0.9329, Adjusted R-squared:  0.9228 
## F-statistic: 92.68 on 3 and 20 DF,  p-value: 6.638e-12
cat("AIC: ", AIC(fit2))
## AIC:  -7.066552
simres <- simulateResiduals(fittedModel = fit2)
plot(simres)

fit3 <- lm(NormCounts2 ~ H2O2*genotype, data = dataset)
print(summary(fit3))
## 
## Call:
## lm(formula = NormCounts2 ~ H2O2 * genotype, data = dataset)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.173733 -0.054959 -0.005033  0.052231  0.190077 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            0.96909    0.05038  19.234 2.27e-14 ***
## H2O2                  -0.49686    0.04275 -11.624 2.38e-10 ***
## genotypeYFP-ALC1       0.06144    0.07125   0.862    0.399    
## H2O2:genotypeYFP-ALC1  0.01404    0.06045   0.232    0.819    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09941 on 20 degrees of freedom
## Multiple R-squared:  0.9301, Adjusted R-squared:  0.9196 
## F-statistic:  88.7 on 3 and 20 DF,  p-value: 9.977e-12
cat("AIC: ", AIC(fit3))
## AIC:  -37.07284
simres <- simulateResiduals(fittedModel = fit3)
plot(simres)

fit4 <- lmer(Counts ~ H2O2*genotype + (1|UID), data = dataset)
print(summary(fit4))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ H2O2 * genotype + (1 | UID)
##    Data: dataset
## 
## REML criterion at convergence: 273.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.35434 -0.61511 -0.05378  0.58692  2.21547 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  UID      (Intercept)  3127     55.92  
##  Residual             31742    178.16  
## Number of obs: 24, groups:  UID, 6
## 
## Fixed effects:
##                       Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)            1330.53      95.89   15.97  13.875 2.51e-10 ***
## H2O2                   -686.47      76.61   16.00  -8.961 1.24e-07 ***
## genotypeYFP-ALC1         60.66     135.61   15.97   0.447    0.661    
## H2O2:genotypeYFP-ALC1    35.28     108.34   16.00   0.326    0.749    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) H2O2   gYFP-A
## H2O2        -0.774              
## gntYFP-ALC1 -0.707  0.547       
## H2O2:YFP-AL  0.547 -0.707 -0.774
cat("AIC: ", AIC(fit4))
## AIC:  285.7402
simres <- simulateResiduals(fittedModel = fit4)
plot(simres)

Quadratic formula

fit5 <- lm(Counts ~ Experiment + poly(H2O2,2)*genotype, data = dataset)
print(summary(fit5))
## 
## Call:
## lm(formula = Counts ~ Experiment + poly(H2O2, 2) * genotype, 
##     data = dataset)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -226.0 -113.0  -48.0  102.0  358.9 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       634.833     77.540   8.187  4.1e-07 ***
## Experimentexp2                     86.125     94.966   0.907    0.378    
## Experimentexp3                      5.875     94.966   0.062    0.951    
## poly(H2O2, 2)1                  -2257.814    268.605  -8.406  2.9e-07 ***
## poly(H2O2, 2)2                    121.603    268.605   0.453    0.657    
## genotypeYFP-ALC1                   94.833     77.540   1.223    0.239    
## poly(H2O2, 2)1:genotypeYFP-ALC1   116.021    379.865   0.305    0.764    
## poly(H2O2, 2)2:genotypeYFP-ALC1  -477.253    379.865  -1.256    0.227    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 189.9 on 16 degrees of freedom
## Multiple R-squared:  0.8966, Adjusted R-squared:  0.8513 
## F-statistic: 19.82 on 7 and 16 DF,  p-value: 9.13e-07
cat("AIC: ", AIC(fit5))
## AIC:  328.218
simres <- simulateResiduals(fittedModel = fit5)
plot(simres)

fit6 <- lm(NormCounts ~ poly(H2O2,2)*genotype, data = dataset)
print(summary(fit6))
## 
## Call:
## lm(formula = NormCounts ~ poly(H2O2, 2) * genotype, data = dataset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.22690 -0.11902 -0.05308  0.10590  0.31224 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      1.000e+00  5.115e-02  19.552 1.42e-13 ***
## poly(H2O2, 2)1                  -3.350e+00  2.506e-01 -13.372 8.67e-11 ***
## poly(H2O2, 2)2                   1.813e-01  2.506e-01   0.724   0.4786    
## genotypeYFP-ALC1                 2.783e-16  7.233e-02   0.000   1.0000    
## poly(H2O2, 2)1:genotypeYFP-ALC1  5.288e-01  3.543e-01   1.492   0.1530    
## poly(H2O2, 2)2:genotypeYFP-ALC1 -6.474e-01  3.543e-01  -1.827   0.0843 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1772 on 18 degrees of freedom
## Multiple R-squared:  0.9451, Adjusted R-squared:  0.9298 
## F-statistic: 61.92 on 5 and 18 DF,  p-value: 1.043e-10
cat("AIC: ", AIC(fit6))
## AIC:  -7.865543
simres <- simulateResiduals(fittedModel = fit6)
plot(simres)

fit7 <- lm(NormCounts2 ~ poly(H2O2,2)*genotype, data = dataset)
print(summary(fit7))
## 
## Call:
## lm(formula = NormCounts2 ~ poly(H2O2, 2) * genotype, data = dataset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.12770 -0.06473 -0.02589  0.05405  0.17573 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      0.48775    0.02716  17.958 6.13e-13 ***
## poly(H2O2, 2)1                  -1.63419    0.13306 -12.282 3.47e-10 ***
## poly(H2O2, 2)2                   0.08843    0.13306   0.665   0.5147    
## genotypeYFP-ALC1                 0.07504    0.03841   1.954   0.0665 .  
## poly(H2O2, 2)1:genotypeYFP-ALC1  0.04617    0.18818   0.245   0.8089    
## poly(H2O2, 2)2:genotypeYFP-ALC1 -0.35076    0.18818  -1.864   0.0787 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09409 on 18 degrees of freedom
## Multiple R-squared:  0.9436, Adjusted R-squared:  0.928 
## F-statistic: 60.28 on 5 and 18 DF,  p-value: 1.308e-10
cat("AIC: ", AIC(fit7))
## AIC:  -38.24457
simres <- simulateResiduals(fittedModel = fit7)
plot(simres)

fit8 <- lmer(Counts ~ poly(H2O2,2)*genotype + (1|UID), data = dataset)
print(summary(fit8))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ poly(H2O2, 2) * genotype + (1 | UID)
##    Data: dataset
## 
## REML criterion at convergence: 241
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2496 -0.5899 -0.2551  0.5390  2.0928 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  UID      (Intercept)  3255     57.05  
##  Residual             31231    176.72  
## Number of obs: 24, groups:  UID, 6
## 
## Fixed effects:
##                                 Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                       665.50      60.73     4.00  10.959 0.000394
## poly(H2O2, 2)1                  -2257.81     249.92    14.00  -9.034 3.24e-07
## poly(H2O2, 2)2                    121.60     249.92    14.00   0.487 0.634098
## genotypeYFP-ALC1                   94.83      85.88     4.00   1.104 0.331429
## poly(H2O2, 2)1:genotypeYFP-ALC1   116.02     353.45    14.00   0.328 0.747575
## poly(H2O2, 2)2:genotypeYFP-ALC1  -477.25     353.45    14.00  -1.350 0.198350
##                                    
## (Intercept)                     ***
## poly(H2O2, 2)1                  ***
## poly(H2O2, 2)2                     
## genotypeYFP-ALC1                   
## poly(H2O2, 2)1:genotypeYFP-ALC1    
## poly(H2O2, 2)2:genotypeYFP-ALC1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(H2O2,2)1 pl(H2O2,2)2 gYFP-A p(H2O2,2)1:
## pl(H2O2,2)1  0.000                                           
## pl(H2O2,2)2  0.000  0.000                                    
## gntYFP-ALC1 -0.707  0.000       0.000                        
## p(H2O2,2)1:  0.000 -0.707       0.000       0.000            
## p(H2O2,2)2:  0.000  0.000      -0.707       0.000  0.000
cat("AIC: ", AIC(fit8))
## AIC:  256.9578
simres <- simulateResiduals(fittedModel = fit8)
plot(simres)

Cubic formula

fit9 <- lm(Counts ~ Experiment + poly(H2O2,3)*genotype, data = dataset)
print(summary(fit9))
## 
## Call:
## lm(formula = Counts ~ Experiment + poly(H2O2, 3) * genotype, 
##     data = dataset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -244.67  -65.79    5.79   39.11  340.21 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       634.833     62.149  10.215 7.18e-08 ***
## Experimentexp2                     86.125     76.117   1.131   0.2769    
## Experimentexp3                      5.875     76.117   0.077   0.9396    
## poly(H2O2, 3)1                  -2257.814    215.292 -10.487 5.17e-08 ***
## poly(H2O2, 3)2                    121.603    215.292   0.565   0.5811    
## poly(H2O2, 3)3                   -410.102    215.292  -1.905   0.0775 .  
## genotypeYFP-ALC1                   94.833     62.149   1.526   0.1493    
## poly(H2O2, 3)1:genotypeYFP-ALC1   116.021    304.468   0.381   0.7089    
## poly(H2O2, 3)2:genotypeYFP-ALC1  -477.253    304.468  -1.567   0.1393    
## poly(H2O2, 3)3:genotypeYFP-ALC1  -170.664    304.468  -0.561   0.5840    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 152.2 on 14 degrees of freedom
## Multiple R-squared:  0.9419, Adjusted R-squared:  0.9045 
## F-statistic:  25.2 on 9 and 14 DF,  p-value: 3.794e-07
cat("AIC: ", AIC(fit9))
## AIC:  318.3933
simres <- simulateResiduals(fittedModel = fit9)
plot(simres)

fit10 <- lm(NormCounts ~ poly(H2O2,3)*genotype, data = dataset)
print(summary(fit10))
## 
## Call:
## lm(formula = NormCounts ~ poly(H2O2, 3) * genotype, data = dataset)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.105261 -0.038084 -0.003288  0.042566  0.114443 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      1.000e+00  2.023e-02  49.429  < 2e-16 ***
## poly(H2O2, 3)1                  -3.350e+00  9.911e-02 -33.805 2.61e-16 ***
## poly(H2O2, 3)2                   1.813e-01  9.911e-02   1.829 0.086062 .  
## poly(H2O2, 3)3                  -6.343e-01  9.911e-02  -6.400 8.79e-06 ***
## genotypeYFP-ALC1                -1.694e-16  2.861e-02   0.000 1.000000    
## poly(H2O2, 3)1:genotypeYFP-ALC1  5.288e-01  1.402e-01   3.773 0.001667 ** 
## poly(H2O2, 3)2:genotypeYFP-ALC1 -6.474e-01  1.402e-01  -4.619 0.000284 ***
## poly(H2O2, 3)3:genotypeYFP-ALC1 -1.210e-01  1.402e-01  -0.863 0.400712    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07008 on 16 degrees of freedom
## Multiple R-squared:  0.9924, Adjusted R-squared:  0.989 
## F-statistic: 296.8 on 7 and 16 DF,  p-value: 1.015e-15
cat("AIC: ", AIC(fit10))
## AIC:  -51.21038
simres <- simulateResiduals(fittedModel = fit10)
plot(simres)

fit11 <- lm(NormCounts2 ~ poly(H2O2,3)*genotype, data = dataset)
print(summary(fit11))
## 
## Call:
## lm(formula = NormCounts2 ~ poly(H2O2, 3) * genotype, data = dataset)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.052134 -0.019338 -0.001694  0.020761  0.055820 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      0.48775    0.01049  46.489  < 2e-16 ***
## poly(H2O2, 3)1                  -1.63419    0.05140 -31.794 6.87e-16 ***
## poly(H2O2, 3)2                   0.08843    0.05140   1.720 0.104632    
## poly(H2O2, 3)3                  -0.30939    0.05140  -6.019 1.79e-05 ***
## genotypeYFP-ALC1                 0.07504    0.01484   5.057 0.000117 ***
## poly(H2O2, 3)1:genotypeYFP-ALC1  0.04617    0.07269   0.635 0.534263    
## poly(H2O2, 3)2:genotypeYFP-ALC1 -0.35076    0.07269  -4.825 0.000186 ***
## poly(H2O2, 3)3:genotypeYFP-ALC1 -0.11570    0.07269  -1.592 0.131009    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03634 on 16 degrees of freedom
## Multiple R-squared:  0.9925, Adjusted R-squared:  0.9893 
## F-statistic: 303.5 on 7 and 16 DF,  p-value: 8.507e-16
cat("AIC: ", AIC(fit11))
## AIC:  -82.72848
simres <- simulateResiduals(fittedModel = fit11)
plot(simres)

fit12 <- lmer(Counts ~ poly(H2O2,3)*genotype + (1|UID), data = dataset)
print(summary(fit12))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ poly(H2O2, 3) * genotype + (1 | UID)
##    Data: dataset
## 
## REML criterion at convergence: 204.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.57978 -0.46992  0.05654  0.30014  2.39500 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  UID      (Intercept)  7219     84.96  
##  Residual             15375    124.00  
## Number of obs: 24, groups:  UID, 6
## 
## Fixed effects:
##                                 Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                       665.50      60.73     4.00  10.959 0.000394
## poly(H2O2, 3)1                  -2257.81     175.36    12.00 -12.876  2.2e-08
## poly(H2O2, 3)2                    121.60     175.36    12.00   0.693 0.501220
## poly(H2O2, 3)3                   -410.10     175.36    12.00  -2.339 0.037475
## genotypeYFP-ALC1                   94.83      85.88     4.00   1.104 0.331430
## poly(H2O2, 3)1:genotypeYFP-ALC1   116.02     247.99    12.00   0.468 0.648277
## poly(H2O2, 3)2:genotypeYFP-ALC1  -477.25     247.99    12.00  -1.924 0.078329
## poly(H2O2, 3)3:genotypeYFP-ALC1  -170.66     247.99    12.00  -0.688 0.504421
##                                    
## (Intercept)                     ***
## poly(H2O2, 3)1                  ***
## poly(H2O2, 3)2                     
## poly(H2O2, 3)3                  *  
## genotypeYFP-ALC1                   
## poly(H2O2, 3)1:genotypeYFP-ALC1    
## poly(H2O2, 3)2:genotypeYFP-ALC1 .  
## poly(H2O2, 3)3:genotypeYFP-ALC1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(H2O2,3)1 pl(H2O2,3)2 pl(H2O2,3)3 gYFP-A p(H2O2,3)1:
## pl(H2O2,3)1  0.000                                                       
## pl(H2O2,3)2  0.000  0.000                                                
## pl(H2O2,3)3  0.000  0.000       0.000                                    
## gntYFP-ALC1 -0.707  0.000       0.000       0.000                        
## p(H2O2,3)1:  0.000 -0.707       0.000       0.000       0.000            
## p(H2O2,3)2:  0.000  0.000      -0.707       0.000       0.000  0.000     
## p(H2O2,3)3:  0.000  0.000       0.000      -0.707       0.000  0.000     
##             p(H2O2,3)2:
## pl(H2O2,3)1            
## pl(H2O2,3)2            
## pl(H2O2,3)3            
## gntYFP-ALC1            
## p(H2O2,3)1:            
## p(H2O2,3)2:            
## p(H2O2,3)3:  0.000
cat("AIC: ", AIC(fit12))
## AIC:  224.6933
simres <- simulateResiduals(fittedModel = fit12)
plot(simres)

Compare Results

ICtab(fit1,fit2,fit3,fit4,
      fit5,fit6,fit7,fit8,
      fit9,fit10,fit11,fit12,
      base=T)
##       AIC   dAIC  df
## fit11 -82.7   0.0 9 
## fit10 -51.2  31.5 9 
## fit7  -38.2  44.5 7 
## fit3  -37.1  45.7 5 
## fit6   -7.9  74.9 7 
## fit2   -7.1  75.7 5 
## fit12 224.7 307.4 10
## fit8  257.0 339.7 8 
## fit4  285.7 368.5 6 
## fit9  318.4 401.1 11
## fit1  327.0 409.7 7 
## fit5  328.2 410.9 9

Final Result

fit <- fit11

output <- coef(summary(fit))
output <- output[grep("H2O2", rownames(output)),]


rownames(output) <- gsub("poly\\(|, [1-3]\\)","", rownames(output) )
rownames(output) <- gsub("genotype",  paste0(" ",levels(dataset$genotype)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$genotype)[1],  sep = " in " )


# suggested result table
kable(output, row.names = T)
Estimate Std. Error t value Pr(>|t|)
H2O21 in WT -1.6341883 0.0513987 -31.7943497 0.0000000
H2O22 in WT 0.0884280 0.0513987 1.7204317 0.1046320
H2O23 in WT -0.3093886 0.0513987 -6.0193863 0.0000179
H2O21: WT vs. YFP-ALC1 0.0461736 0.0726887 0.6352232 0.5342627
H2O22: WT vs. YFP-ALC1 -0.3507551 0.0726887 -4.8254386 0.0001864
H2O23: WT vs. YFP-ALC1 -0.1157002 0.0726887 -1.5917212 0.1310090
write.table(output, file = "FigureS5N_Stats.txt", quote = F, sep = "\t", row.names = T, col.names = NA)

Anova

fit11a <- lm(NormCounts2 ~ poly(H2O2,3)*genotype, data = dataset)
fit11b <- lm(NormCounts2 ~ poly(H2O2,3)+genotype, data = dataset)

# anova table
anova(fit11a, fit11b)
## Analysis of Variance Table
## 
## Model 1: NormCounts2 ~ poly(H2O2, 3) * genotype
## Model 2: NormCounts2 ~ poly(H2O2, 3) + genotype
##   Res.Df      RSS Df Sum of Sq      F   Pr(>F)   
## 1     16 0.021135                                
## 2     19 0.055772 -3 -0.034637 8.7406 0.001156 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1